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, _ Abstract 

We show how the path-integral formulation of quantum statistical mechanics can be used to construct practical ab initio 
techniques for computing the chemical potential of molecules adsorbed on surfaces, with full inclusion of quantum nuclear 
" ^ effects. The techniques we describe are based on the computation of the potential of mean force on a chosen molecule, 
and generalise the techniques developed recently for classical nuclei. We present practical calculations based on density 
' I functional theory with a generalised-gradient exchange-correlation functional for the case of H2O on the MgO (001) surface 
^ at low coverage. We note that the very high vibrational frequencies of the H2O molecule would normally require very large 
numbers of time slices (beads) in path-integral calculations, but we show that this requirement can be dramatically reduced 
• by employing the idea of thermodynamic integration with respect to the number of beads. The validity and correctness 
d of our path-integral calculations on the H20/MgO (001) system are demonstrated by supporting calculations on a set of 
simple model systems for which quantum contributions to the free energy are known exactly from analytic arguments. 

I 

C 1 Introduction 
O 

I ^ i ln a recent paper referred to here as paper I, we described an ab initio scheme for calculating the chemical potential 
and other thermodynamic properties of systems of adsorbed molecules, with full inclusion of entropy effects. We presented 
illustrative calculations based on density functional theory (DFT) for H2O at low coverage on the MgO (001) surface, 
and showed how the chemical potential, and hence the desorption rate can be calculated without statistical-mechanical 
, approximations, apart from the treatment of the nuclei as classical particles. That initial study was intended as a first 
step in developing ab initio methods for calculating the thermodynamic properties of adsorbate systems with better than 
chemical accuracy, i.e. better than 1 kcal/mol. As a contribution to these developments, we present here a practical scheme 
for including quantum nuclear effects in the ab initio statistical mechanics of adsorbates, and we report here some results 
for H2O on MgO (001) at low coverage. 

Quantum effects are likely to have a very significant effect on the properties of some adsorbates, particularly those 
containing protons, such as H2O. The zero-point energies ^Tiw of the symmetric and antisymmetric stretching modes and 
the bending mode of the H2O molecules are 227, 233 and 99 meV respectively (5.23, 5.37 and 2.28 kcal/mol) [2l. Since 
partial or complete dissociation of adsorbed H2O occurs quite commonly [SI lU [5l |6l [7] , it is clear that inclusion of quantum 

r \ nuclear effects will sometimes be essential for accurate calculations. Needless to say, the future achievement of chemical 
accuracy will also depend crucially on the development of ab initio methods which go beyond conventional DFT, and that 
is currently the object of intensive research [51 1^ [TUl 111] . However, the availability of improved ab initio methods is not 
our immediate concern here. 

We emphasised in paper I that we aim to avoid statistical-mechanical approximations, and in particular we do not allow 
ourselves to resort to harmonic approximations, since these will not be satisfactory for the disordered adsorbates that are 
of ultimate interest. The overall strategy we shall present here closely resembles that of paper I, the main new feature 
being that the path-integral formulation of quantum mechanics is used to treat the nuclei. We noted in paper I that the 
chemical potential determines the ratio between the surface adsorbate density and the gas-phase density with which it is in 
equilibrium, and that this ratio can be expressed as an integral over the normalised density distribution y{z) of molecules 
as a function of height z above the surface. Since y{z) can be expressed in terms of a potential of mean force (t>{z), it can 
be computed from a sequence of thermal-equilibrium simulations in each of which the mean force on a chosen molecules is 
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calculated, with the height z of this molecule above the surface being constrained to have a fixed value. As we shall see, 
these relationships remain true in path-integral theory. 

Path-integral statistical mechanics is based on the well known isomorphism between a quantum system of N particles 
and a classical system of N cyclic chains, each consisting of P beads (time slices), with neighbouring beads in each chain 
being coupled by harmonic springs [T5( [T21 EH UHl [IHl HZl HH] ■ The meaning of this is that each cyclic chain represents 
the propagation of an individual quantum particle in imaginary time r from r = to r = fih. The ratio between the 
densities of adsorbed and gas-phase molecules, and the expression for this ratio in terms of a probability distribution y{^z) 
retain their validity when the nuclei are described by quantum mechanics, and we shall show that y{z) can be expressed in 
terms of a potential of mean force 0(z) which can be derived from appropriately constrained simulations on the system of 
cyclic chains given by the path-integral isomorphism. In this sense, the generalisation of the methods of paper I to include 
quantum nuclear effects is in principle straightforward. 

We shall use the example of H2O on MgO (001) to show that it is practically feasible to calculate the chemical potential 
of an adsorbate using ah initio path-integral methods. The illustrative calculations are for the case of low coverage, i.e. 
isolated molecules, because we are not yet able to perform the statistical sampling at higher coverages. For isolated H2O 
molecules on MgO (001), the changes of vibrational frequencies on adsorption are not large, so that quantum nuclear 
effects do not lead to great changes of the chemical potential. Nevertheless, the technical challenges of the calculations are 
interesting and substantial, because the high vibrational frequencies mean that the number of path-integral beads needs 
to be very high to achieve the required accuracy. To meet this challenge, we had to develop a new technique, in which 
thermodynamic integration is performed with respect to the number of beads. A significant part of the paper is concerned 
with the tests we have developed to demonstrate that our path-integral techniques work correctly. 

The rest of the paper is organised as follows. The techniques themselves are laid out in Sec. [2j we start with a 
brief reminder of the methods developed in paper I, and a summary of basic path-integral theory; then we explain how 
this theory can be used to create a path-integral generalisation of paper I that correctly yields the chemical potential 
with full inclusion of quantum nuclear effects; techniques for improving the computational efficiency and for performing 
thermodynamic integration with respect to the number of beads are also outlined; at the end of Sec. [2] we note the DFT 
techniques used in our practical calculations. Sec. [3] describes the tests we have performed on the path-integral machinery; 
these tests consist of simulations on a set of harmonic vibrational models involving the H2O molecule acted on by various 
external fields, for which the results can be checked against values known from analytical formulas. The calculations on 
the chemical potential of the H2O molecule adsorbed at finite temperature on MgO (001) are presented in SecJ4] Finally, 
the implications of the work and the prospects for future applications and developments are discussed in Sec. [5] 



2 Techniques 

2.1 Chemical potential of adsorbed molecules: classical approximation 

The classical theory of paper I is expressed in terms of the spatial distribution p{z) of molecules when there is full thermal 
equilibrium between the gas phase and the adsorbate. This distribution is defined so that p{z) dr is the probability of 
finding a molecule in volume element dr at height z above the surface, averaged over positions in the plane parallel to the 
surface. The detailed form of p(z) depends on which point within the molecule (e.g. centre of mass, position of a specified 
atom in the molecule, etc.) is chosen to construct the distribution; we refer to this as the "monitor point" [1]. In practice, 
it is convenient to work with the normalised function y{z) = p(z)/po, where po is the mean density of molecules in the gas 
phase. The ratio of the density a of adsorbate molecules (number per unit area on the surface) and the gas-phase density 
Po is: 

dzy(z) , (1) 

-00 

where zq is the height above the surface below which molecules are counted as being "adsorbed" . This ratio is also related 
to the excess chemical potentials of gas-phase and adsorbed molecules: cr/pa — dexp[/3A/j,'''(f7, T)], where Ap^ = M|as~A'Ids 
is the difference of excess chemical potentials, defined in paper 1, (3 = l/k^T and d is a length chosen for the purpose 
of defining /i^ds- We refer to A/i^ as the "excess chemical potential difference" (ECPD). In principle, the choice of zq 
affects cr, and hence Ap^; but in most situations of interest y{z) for z < zg is enormously greater than the gas-phase value 
y{z — ^ cxi) — 1, so that in practice the dependence of cr/po on zq can be neglected. For the same reason, the dependence 
of cr/po and A/i^ on the choice of monitor point is extremely weak. 

Since y{z) » 1 in the adsorption region z < zq, the computation of y{z) is a "rare event" problem. To overcome 
this problem, y{z) is expressed as y(z) = exp(— /30(z)), where the potential of mean force (j){z) represents the free energy 
of the entire system when the height of the molecule above the surface is constrained to have value z, with the condition 
(j){z 00) — 0. To compute (piz), we use the fact that (-Fz)z = —d4>/dz is the z-component of the thermal average force 
acting on the molecule when its height is constrained to be z. In the practical calculations of paper I, ab initio molecular 
dynamics (m.d.) simulation was used to perform the constrained statistical sampling in the canonical ensemble, and hence 
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the computation of {Fz)z- Integration of {Fz)z then yields (t){z) and hence y{z). Related techniques have been developed 
by other groups (see e.g. Ref. [I9]). 



2.2 Basic path-integral methods 

We summarise briefly the path-integral background needed in this work; for details, the many reviews should be con- 
sulted [m [131 mi [ini [m [TS] . The Helmholtz free energy F of a general system at temperature T is F = —k^TlnZ, where 
the partition function is Z = Tr [exp(— with j3 = l/k-oT and H — T + V the Hamiltonian of the system. Here T is 
the operator for the sum of kinetic energies of all the nuclei, and V = [/({r^}) in the present work represents the electronic 
ground-state energy of the system when the nuclei are fixed at positions r.^ (z = 1. 2, . . . N). 

As usual in path-integral theory, Z is rewritten as Z = Tr [(/5(/3/P))^] , where p{f3/P) = exp{—{f3/ P)H) is the 
propagator for imaginary time f3h/P, which is approximated by the Trotter short-time formula: 

p(/3/P)^e-^^/2^po(/3/P)e"^^/^^, (2) 
with pq{P / P) the free-particle propagator, given in coordinate representation by: 



Po(R,R';/3/P) = (R|po(/3/P)|R') = Apexp 



miP 



(3) 



Here, R and R' are points in the configuration space of the whole system, specified by the collections of positions {r^} and 
{r^}, and the prefactor Ap is given by: 

with rrii the mass of nucleus i. The resulting approximation to the partition function, denoted here by Zp, is then: 

„P-1 N 

Zp - {Apf / n n exp[-/3(rp({r, J) -I- Vp{{v,,})] , (5) 

s=0 i=l 

where r^s is the position of nucleus i at time-slice s, and Tp and Vp are given by: 

p-i N ^ 

Tp{{r,s}) = ^ ^ -Kpj|r,,+i - r,s|^ 

P-1 

P 



Vp{{v.s}) = \,Y.U{{v^s}) ■ (6) 



s=0 



The spring constants np^i are given by Kp^i = mi P/P'^Ti'. The formula for Zp expresses the widely used approximation to 
the partition function of the Ai"-particle quantum system in terms of the partition function of a classical system of cyclic 
chains, each consisting of P "beads" (time slices), with neighbouring beads in each cyclic chain being coupled by harmonic 
springs. The exact partition function is recovered in the limit P — ?> cxd. 

The spatial distributions of the positions of beads belonging to a particular time-slice in the classical system of cyclic 
chains have a special significance, because they are equal (in the P — ^ cxd limit) to the corresponding equal-time spatial 
distributions in the quantum system. For example, the probability density of finding a particular nucleus at a particular 
spatial point r in the quantum system is equal to the probability density of finding a chosen bead on the cyclic chain 
representing that nucleus at point r. This means that appropriate potentials of mean force in the classical cyclic-chain 
system allow us to calculate spatial probability distributions in the quantum system. 



2.3 Constraints and mean force in path-integral scheme 

In the case treated later of H2O adsorbed on MgO (001), we could choose to take as monitor point the position of the 
O atom in the H2O molecule. In that case, the normalised height distribution y{z) of this monitor point in the quantum 
system is equal to the height distribution y{z) of a chosen bead on the cyclic chain representing that O atom in the classical 
cyclic-chain system. Alternatively, if we took as monitor point the centre of mass of the H2O molecule, i.e. the position 
j,cm _ (^jjio^^i + TnY{V2 + mur^) / {mo + 2mH) in an obvious notation, then the height distribution y{z) of the z-component 
z"^ of this position in the quantum system is equal to the height distribution y{z) of the quantity z*^™ constructed from the 
O and H positions, all at the same chosen time-slice in the classical cyclic-chain system. Just as in the classical theory of 
paper I, the detailed form of y{z) will depend on the choice of monitor point, but the integral JJ'^ dzy(z) will not depend 
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appreciably on this choice (see above). In the quantum case, there are also other possible ways of constructing monitor 
points. One possibility is that we monitor the centre of mass of the centroids of the cyclic paths of the O and the two H 
atoms of the H2O molecule. (The centroid of the cyclic chain of nucleus i is the average = SfLi ''is-) If we do this, 
then y{z) for the classical system of cyclic chains does not necessarily correspond to any obervable of the quantum system. 
Nevertheless, the integral Jf^ dzy{z) will still be the same as for other choices of monitor point. This choice of "centroid 
centre of mass" is the one that we use in the present work. 

For any of the foregoing choices of monitor point, the distribution y{z) can be computed from the potential of mean 
force 4'{z), defined in the usual way as y{z) — exp(— /30(z)), and the derivative —d(j)/dz is {Fz)z, the mean value of the 
z-component of the force on the monitor point, when the z-component of the position of this monitor point is constrained 
to have a particular value. 



2.4 Thermal sampling 

When performing constrained thermal-equilibrium simulations of the cyclic chain system in order to calculate quantities 
such as {Fz)z, we can choose to use any simulation method that samples configuration space in the required manner. As in 
paper I, we use molecular dynamics (m.d.) simulation with Nose [50] and Andersen [5T] thermostats to perform sampling 
according to the canonical ensemble. We noted in paper I that with this approach we are free to choose the m.d. masses 
of the particles in any way that helps to improve the sampling efficiency, since spatial averages do not depend on these 
masses. When we use m.d. to perform thermal sampling of the system of cyclic chains, it is therefore essential to distinguish 
between the real nuclear masses rrii appearing in the spring constants Hi from the m.d. "sampling masses" employed to 
sample the configuration space. In paper I, the sampling masses were chosen as Mh = 8, Mq = 16 and MMg = 24 a.u. 

It is well known that if m.d. simulation is used to perform thermal sampling in the path-integral technique, with each 
bead on every chain being give a chosen mass, then the sampling efficiency is very poor, because of the very wide frequency 
spread of the vibration modes of each chain. The method we use to overcome this problem, which resembles the techniques 
used by other workers (see e.g. Ref. |22j). employs a mass matrix on each cyclic chain, so that the Hamiltonian governing 
the m.d. evolution is: 

^ N P-lP-l ^ N P-l 

^=2^^^^'"'^*'*'^'* ^ 2 51X1 '^^I'^'^+i 

i=l s=0 t=0 i=l s=0 

P-l 

s=Q 

where A^st are the elements of the (positive definite) inverse mass matrix for the cyclic chain representing nucleus i. Each 
mass matrix is chosen so that the vibration frequencies of all modes of each chain are similar, and so that the total mass 
of each chain has a value that would be suitable for the sampling mass in the classical system. This is achieved if we take: 

{CjP)cos[2nn{s-t)/P] 
^„ iQ/P) + 2k,{1 - cos(2^n/P)) ' ^ ^ 

The total sampling mass Mi of each chain is chosen exactly as in the classical simulations of paper I (see above). The 
spring constants Ci are chosen so that (Ci/Mi)^^^ is a typical frequency characterising the dynamics of nucleus i in the 
classical simulations performed with masses M^. 



2.5 Switching number of beads 

In some of our calculations, it is helpful to obtain the free energy by successive approximations. First, we calculate 
Fp = —k^ThiZp with a number of beads P which is known to be insufficient. We then obtain a more accurate value F2P 
by adding the difference F2P ~ Fp. If the accuracy is still insufficient, we repeat the process. There is a simple way of 
using thermodynamic integration to calculate the difference F2P — Fp, which we now explain. 

We define a generalised form of the free energy F2P for the system having 2P time slices by introducing weights Wg 
(s = 0, . . . 2P — 1) into the definition of the potential energy V of cqn ([6]): 

2P-1 



^2p({r,.}) - ^ E ^sUi{r,,}) 



2P 



s=0 



The partition function Z2P and the free energy F2P are now functions of the Wg- We choose the 
Wg = 1 + X for even s and Wg = 1 — A for odd s, so that the partition function Z2P is given by: 



(9) 



to have the form 



Z2P — {A2p) 



2P 



2P-1 N 



n n^r. 



exp 



( T2p{{r,g}) + ^ E (1 + i-iyX)Ui{r,g}) 



2P-1 



2P 



(10) 
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Clearly, Z2P and hence the free energy Fip depend on A. From eqn (10), we readily obtam a formula for the variation of 
Fip with A: 

aF2p/aA = -(^Xi(f^({'-'2s+i})-^({r.2.})) . (11) 

Now the quantity F2p{X — 0) is just the usual approximation to the free energy calculated with 2P beads, as is clear 
from eqn (10). On the other hand, F2p{X = 1) is identical to Fp{X — 0), which is the normal free energy calculated with 
P beads. To see this, we note from eqn ([To| that for A = 1 the potential term in the Boltzmann factor becomes: 



2P-1 



P-1 



2^ 5^ (1 + (-l)^)C/({r,,}) ^^Y. f^({r'2.}) , 



(12) 



which is the potential term that appears in Fp(X = 0). But this means that we can integrate over all the positions r^g for 
odd s, since these no longer appear in the potential term. In doing the integrations, we use the fact that: 



2 f ^ r ^ 1 

{A2pf / J|dri2s+i exp -/3^-K2P, 

1=1 L i=l 

^ 1 



Fi2s+i - ri2s| exp 



^ 1 



K2P,i|ri2s+2 — ri2s + l| 



i=l 



Ap exp 



I^PA\'ri2s+2 — Yi2s 



i=l 



(13) 



which expresses the fact that /5o(/3/(2P))/5o(/?/(2P)) — /5o(/3/P). On using eqns (12) and (13) in eqn (10), we obtain 
F2P(A = 1) = Fp{X = 0). 



It follows from this, and from eqn (fill), ^^^^ 

»i 



F2P -PP = ly^i^ 5]V({r.2s+i}) - U{{v.,2s}) 



(14) 



/ A 

This is the formula we use to calculate accurate values of free energy by successive approximation. 



2.6 Ah initio techniques 

Our practical simulations, like those of paper I, were made with the projector-augmented-wave implementation of DFT [331 
[231, using the VASP code [25,. The DFT computations on all the P path-integral images of the system (i.e. beads, or time- 
slices) are performed in parallel, with the operations for each image also distributed over groups of processors (typically, 
there are 16 processors allocated to each image for our full calculations on the MgO slab with or without II2O). As in 
paper I, we use the PEE exchange-correlation functional ^26j, with a plane- wave cut-off of 400 eV and an augmentation- 
charge cut-off of 605 eV. 



3 Testing the path-integral machinery 

To verify that the path-integral techniques work correctly, we have used them to calculate the free energy changes when 
various external potentials are applied to an isolated II2O molecule in free space. We have designed these potentials so 
that they mimic the changes that occur when a molecule passes from free space to the adsorbed state: the formation of 
vibrational and librational modes from free translations and rotations, and the alteration of intramolecular vibration modes. 
The potentials are designed so that the changes of free energy are known exactly. These tests also give us information 
about the number of beads P needed to obtain results of specified accuracy. All the tests are performed with appropriately 



modified versions of the VASP code, with the technical settings mentioned above (Sec. 2.6) 



3.1 Vibration of centre of mass 

To mimic the conversion of free translations to vibrations, we introduce a potential that acts only the centre of mass of the 
molecule. This potential leaves molecular rotations and internal vibrations completely unchanged. We use the notation 
for the positions of the three atoms in the molecule, with « = 1, 2, 3 corresponding to the O atom and the two H atoms 
respectively. The position of the physical centre of mass is r™^ = (mofi + TOH(r2 + r3))/M, where M = mo + 2mH is the 
total mass of the molecule. The total energy J7o(ri, r2, rs) of the molecule in free space depends only on relative positions, 
such as Y2 — Vi and ra — ri, since it is translationally invariant. We add to C/q an external potential C/cxt(2'^'") depending 
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only on the z-coordinate of the centre of mass. We want Ucxt{z™) to consist of a harmonic well of spring constant a and 
depth Vq, and we want it to go to zero as z — > oo. To obtain this, we adopt the form: 

Uext{z) = -Vo + ^a{z ~ Zq)'^ , z < Zi 

= -^f3{z~Z2)^, Zl < Z < Z2 

= , Z2<z. (15) 

The constants are chosen so that Ucxt and its first derivative are continuous at 2 = zi. In this model, the rotational and 
internal- vibrational properties of the molecule are completely unchanged, and the only effect of Ucxt is to produce harmonic 
oscillations of the centre of mass in the z-direction when z < zi. We take the potential parameters and T to be such that, 
when the molecule is in thermal equilibrium in the potential well, the probability of finding the centre of mass outside the 
region z < zi is negligible. Then the values of A/i^ in the quantum and classical cases are exactly A:BTln[2sinh(?ia;/2fcBr)] 
and k^T ln{tiLL) / k^T) , so that their difference (quantum minus classical) is kBThi[smh{huj/2kBT)/{hLL)/2k-BT)], where 
uj = yfajM is the oscillation frequency of the centre of mass. 

In the classical case, simulations are superfluous, since if the z-component of the centre of mass position is constrained 
to be z, then the force on the centre of mass is simply Fz — —dUcyit/ dz. In the quantum case, the coordinate we choose to 
constrain is the z-component of the centre of mass of the centroids z"^™ = mozi 4-mH(z2 -l-Z3))/M. Now it is readily shown 
that for a perfectly harmonic system in thermal equilibrium, the probability distribution of the path ccntroid is identical 
to that of the corresponding classical system. This implies that the mean force is also identical, and that the potentials of 
mean force in the quantum and classical cases can differ by at most a constant. In the present case of a piecewise harmonic 
[/ext, the mean force is therefore the same in the quantum and classical systems, except for those z'^™ values for which 
beads are distributed on both sides of the "join" position zi. Contributions to the quantum-classical difference of A/i^^ 
therefore come only from this join region. 

We have performed practical tests on this model for a variety of parameter values. An example is the case a = 
31.975 eV/A2, which gives a centre-of-mass oscillation frequency of 20.8 THz. With |zi — zq I = 0.5 A (we take {3 = a/25), 
the centre of mass is almost completely confined to the region z < zi in both the quantum and classical cases at T = 100 K. 
In this example, the exact quantum-classical difference of A/i^^ is 23.2 meV. Our path-integral calculations of mean force 
for this case, performed with P = 8 (spot checks with P = 16 showed no appreciable differences), gave a difference of 
23.2 ± 0.5 meV, in perfect agreement with the exact value. 



3.2 Libration of normal to molecular plane 

To illustrate the conversion of free rotations to librations, we take a potential that acts on the vector normal to the 
molecular plane of the II2O molecule. A suitable vector is: 

Q± = (r2-ri) X (r3-ri) , (16) 

and we define the potential as: 

V,xt = \a[nl+nl) , (17) 

where the unit vector n is Q_l/|Q_l|. For a = 0, the molecule rotates freely, but for large positive a the normal executes small 
oscillations about the z-axis. If a is large enough, these oscillations are harmonic, and their frequencies are lui — (a//i)^/^ 
and UJ2 — (a/^2)^^^j where Ii and I2 are the moments of inertia for rotations about the two symmetry axes in the molecular 
plane. In-plane rotations (moment of inertia /a) are unaffected. The free energy in the presence of T^xt is then the free 
energy of the two librations plus that of the free in-plane rotation. For the free energy in the absence of Vcxt, we use the 
standard expression for the free energy of an asymmetric top [27) . The free energy increase AP caused by switching on 
T4xt is thus known almost exactly. For our tests, we take a = 2.00 eV and T = 100 K; under these conditions, the two 
librational modes are almost exactly in the quantum ground state. The free energy increase is AP = 111 meV. 

We test the path-integral methods using a series of simulations with scaled external potential AT4xt, with A going 
from to 1. The thermal average of V^xt in the path-integral system was calculated using A values of 0.0, 0.03, 0.06, 
0.125, 0.25, 0.5 1.0, the reason for taking closely spaced A values near A = being that the thermal average varies 
rapidly with A in this region (see Table [I]). Numerical integration gives AP = 102 ± 1 meV, and AP = 106 ± 1 meV 
using P = 8 and P = 16 respectively. The small free energy difference obtained between P = 8 and P = 16 was also 



reproduced by directly calculating it using eqn (14), which gives APig — APg = 5±1 meV. Using eqn (14), we also obtained 
AP32 — AP16 = 1.5 ±0.5 meV, so that our best estimate for AP is 108.5 ±2 meV, which agrees with the exact value within 
the statistical errors. 
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P 


A 


8 


16 


0.0 


680 ± 11 


664 ± 12 


0.03 


405 ± 12 


443 ± 15 


0.06 


256 ± 8 


247 ± 10 


0.125 


153 ±4 


144 ±5 


0.25 


97 ± 2 


101 ± 3 


0.5 


65 ± 1 


70 ± 2 


1.0 


41 ± 1 


47 ± 1 


AF 


102 ± 1 


106 ± 1 



Table 1: Thermal average {Vext)x as function of A and free energy change AF. Units are meV. 



p 


AF(meV) 


96 


78 ± 1.5 


64 


75 ±1.5 


32 


59 ±1.5 


16 


41 ± 1 


8 


21 ± 1.5 



Table 2: Free energy of the water molecule (in meV) in the external potential Zcxt at 100 K calculated using path integral 
simulations with different values of P. 



3.3 Alteration of asymmetric stretch frequency 

Surface adsorption will alter the intramolecular vibration frequencies, and to mimic this effect we apply an external potential 
having the form: 

Zo.t - , (18) 

where ^ = (r2 +r3 — 2ri) • (ra — r2). Clearly, this has no effect on translations or rotations. By symmetry, ^ remains zero in 
the presence of symmetric stretch and bond-angle oscillations. The only mode that is affected by Zcxt is thus the asymmetric 
stretch mode. In the absence of Z^xt, our DFT calculations give the frequency aJa/27r — 114.8 THz. We have done tests 
with the value n — 5.0236 eV/A^, which gives the increased frequency uji^/2n — 151.2 THz. At T = 100 K, this vibrational 
mode is in its quantum ground state, so that the free energy change caused by Zcxt is AF — ^^{uja ~ i^a) — 77.3 meV. 

In our path-integral simulations, we calculate the thermal average of Zcxt as a function of A in the presence of the 
scaled potential XZcxt, with A going from to 1. These calculations were done with several different numbers of time-slices 
P, going from 8 to 96. The computed results for the free energy change AF (Table [2| show that P = 96 gives essentially 
exact results, and P = 64 is also acceptable, but smaller values give seriously inaccurate results. This indicates that if we 
wished to calculate the chemical potential of adsorbed H2O by brute force path-integral simulation, we would have to use 
at least P = 64, which would be exceedingly expensive. 

This gives us an opportunity to test our technique of thermodynamic integration with respect to P. To apply this, we 
first calculate the free energy change as A goes from to 1, using P = 8. The resulting value of AF is in error by nearly 
50 meV (Table [2]). To correct this, we now calculate the free energy changes P(P) — P(P/2), both for the free molecule 
and for the molecule acted on by Zcxt- The resulting correction gives us AP(8) -I- AP(8 — > 64) = 70 ± 6 meV, which agrees 
within statistical error with the directly calculated value AP(64) = 75 ± 1.5 meV. 

3.4 Other tests 

As a further test of librational effects, we have compared path-integral calculations with quasi-exact results for the case 
where potentials act both on the normal to the molecular plane and on the direction of the molecular bisector. This is 
achieved by adding to the potential Vcxt of Sec. |3.2| a potential Wcxt defined by: 

Wcxt = \iql , (19) 

where is the x-component of the unit vector given by q = (r2 -f — 2ri)/|r2 -I- r3 — 2ri|. The tests were done with 
7 = 2.00 eV, which gives a frequency of ^ 16 THz for in-plane librations. The potential Wext also induces a small change 
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Figure 1: Mean force {Fz)z and potential of mean force (t){z) of H2O on MgO (001) calculated here using path-integral 
simulation with P = 8 beads (diamonds), compared with the classical results of paper I (circles). Units are A and eV. 
Statistical error bars are smaller than the size of the points. 



in the frequency of the asymmetric stretch mode. As before, the free energy change caused by switching on V^oxt + W^oxt 
agrees within statistical error with the quasi-exact value. 

In addition to all these tests on the ah initio free molecule, we have also made a completely different kind of test 
involving PMF calculations on the molecule as it is brought onto the surface of the rigid MgO slab (see Appendix). The 
methods are the same as those used for the full path-integral calculations presented in the next Section. 



4 Full path-integral calculations for H2O on MgO (001) 

The PMF techniques of Sec. [2] are now applied to the full calculation of the ECPD A/i^ = //^.g^g — Acids' both with and 
without nuclear quantum effects. We denote the difference of A/i^ values (quantum minus classical) by (A/i^)qc. We 
know from the tests of Sec. [s] that in order to obtain values for (A/i^^)qc that are correct to better than ^ 5 meV, the 
number of beads P must be no less than 64. Our strategy is to perform all PMF calculations with P = 8, and then to 
correct the results by performing thermodynamic integration with respect to the number of beads, for the free molecule, 
for the slab-molecule system, and for the clean slab. All the parameters of the simulated system were set exactly as in 
paper I: we use a 3-layer 2x2 slab (16 ions per layer, for a total of 48 ions per repeating unit of the MgO slab), with a 
vacuum gap of 12.7 A, a lattice parameter of 4.23 A (the PBE bulk equilibrium value at T = K), and F-point electronic 
sampling. We use a time-step of 1 fs, and m.d. sampling masses of 24.3, 16.0 and 8.0 a.u. for Mg, O and H, in conjunction 



with the mass-matrix technique (Sec. 2.4). With these parameters, the conservation of the constant of motion was of 
the same quality as that obtained in the classical simulation of paper I. The calculations were performed at T = 100 K, 
with an Andersen thermostat |21| . and constrained simulations were made at 13 different values of height z. We find that 
simulations of 20 ps suffice to give a statistical error on the PMF at its minimum of less than 6 meV, which, as shown in 
paper I (Appendix A), is approximately equal to the error in the ECPD A/i^. 

The mean forces {Fz)z and the PMFs 0(2;) calculated with P = 8 and P = 1 (classical case) are shown in Fig. [l] We 
see that quantum effects are small, but are easily detectable within the statistical errors. In particular, we notice that the 
classical mean force is slightly lower than the quantum value almost everywhere, apart from the region 4.8 < z < 5.2 A, 
where it is higher. By integrating the PMFs, we obtain a quantum-classical difference of EPCD of (A/i^^)qc = —22 ±7 meV. 

In computing the free energy changes on going from P = 8 to P = 64 for the isolated molecule and slab-molecule system, 
we used 4, 3 and 2 values of A for the 8 — > 16, 16 — > 32 and 32 — >■ 64 contributions (simulations at A = are unnecessary, 
since {Up — U2p)\=o — 0). For the clean slab, we needed only two values of A, since the P = 8 — 64 corrections for 
this system are rather small. In Fig. [2] we show the values of {Up — U2p)\ as a function of A for the three systems. We 
see that most of the nuclear quantum contribution comes from the step 8 16, but the other two contributions are not 
negligible. However, when we combine the contributions from the three systems, the three contributions 8 — >■ 16, 16 — >■ 32 
and 32 — >• 64 are all very small, and the total change of A/i^ (64 beads minus 8 beads) amounts to only — 8 ± 3 meV. 
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Figure 2: Computed results for the quantity {Up — U2p)x (eV units) used to calculate the difference Fp — F2P between 
path- integral free energies for numbers of beads P and 2P. Panels (a), (b) and (c) show results for the slab- molecule system, 
the free molecule, and the bare slab respectively. Results for P = 8, 16 and 32 are shown by squares, circles and diamonds. 
Statistical error bars are smaller than the size of the points. 



Including now the correction due to going from P = 8 to 64 beads, the final result for the nuclear quantum correction 
to the ECPD comes to (A/xt)^^ = -30 ± 8 meV. 



5 Discussion and conclusions 

The main purpose of this work has been to establish the practical possibihty of calculating the chemical potential of 
adsorbed molecules using ab initio methods with full inclusion of nuclear quantum effects, but without resorting to harmonic 
approximations. We have shown that within the path-integral formulation of quantum mechanics this can be achieved by 
a generalisation of the ab initio methods that have already been used for the case of classical nuclei. In particular, the 
method based on computing the mean force acting on a single molecule in a series of constrained simulations generalises 
naturally to the path-integral framework. Our practical calculations on the chemical potential of H2O on MgO (001) at 
low coverage show that available computer power suffices to reduce statistical errors to ~ 10 meV, if required. However, 
this is achievable only thanks to a techniques we have introduced for performing thermodynamic integration with respect 
to the number of beads (time-slices) . 

Naturally, future calculations of this kind will be of greatest interest when the nuclear quantum effects are large, and 
this will typically be true when very high intra-molecular vibrational frequencies are strongly affected by surface adsorption. 
Many examples of this are known, among which are water adsorption on both oxide and metal surfaces. For example, it has 
recently been shown that for water layers on some transition-metal surfaces the traditional distinction between covalent and 
hydrogen bonds can be partially or almost entirely lost [35]. However, for the H20/MgO (001) case studied here, nuclear 
quantum effects turn out to have a rather small net effect on the chemical potential, changing it by only 30 meV. This is 
not because the the zero-point energies are small (they are not), but because the intra-molecular vibrational frequencies of 
H2O are not greatly changed on adsorption, and the effects of the changes are partly cancelled by the zero-point energies 
of the new molecule-surface modes created from translational and rotational modes of the free molecule. However, this will 
not be the case for many other water-adsorption systems, and we believe there is ample scope for interesting applications 
of our methods. 

The practical calculations of this paper and paper I were on the H20/MgO (001) system at very low coverage, where 
interactions between adsorbed molecules can be ignored. However, we focused on this low-density limit out of necessity 
rather than choice. In most practical situations, intermolecular interactions play a major role, and certainly have a strong 
effect on thermal desorption rates and on the surface coverage in thermal equilibrium with a given gas-phase pressure. 
For ab initio statistical mechanics, calculations away from the low-density limit present two kinds of problem. First, the 
simulated systems need to be larger, in order to reduce system-size errors to an acceptable level. Second, memory times, 
and hence sampling times usually become much longer, because molecules tend to hinder each other's translational and 
rotational motion. We know from our own work with empirical interaction models for H20/MgO (001) [33ll34] that, if we 
use the PMF approach, simulation times need to be at least 10 times as long at higher coverage than in the low-density limit. 
Nevertheless, with empirical interaction models, the calculation of adsorption isotherms (coverage as a function of chemical 
potential at fixed T) has been routinely practiced for many years, using grand canonical Monte Carlo (GCMC) [3SJI37] and 
other techniques. It is clearly an important challenge for the immediate future to develop the techniques that will allow 
the same to be done with ab initio methods, and ultimately with ab initio path-integral methods. This will not necessarily 
be done with GCMC, though, because there may be more effective strategies. 
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Before concluding, we return briefly to the issue of improving ab initio accuracy, which we noted in the Introduction. 
Improvements are much needed, because DFT energetics of surface systems often depends significantly on the approximation 
used for exchange-correlation energy. Increasingly effective methods for applying highly-correlated quantum chemistry to 
extended systems [Ml [Ml H UHl IMl IS] , the growing power of quantum Monte Carlo techniques [55 1 [M l liO l li H [51 H5] . 
and progress with improved van der Waals DFT functionals [44J give promising signs that ab initio calculations on surface 
adsorbates with chemical accuracy or better are now coming within reach. These improvements will enhance the importance 
of being able to account for nuclear quantum corrections. 

In conclusion, we have presented ab initio path-integral techniques which allow the calculation of the chemical potential 
of adsorbate molecules in thermal equilibrium, with the nuclei treated fully quantum mechanically. We have verified that 
the techniques give correct results by applying them to a set of models where the quantum corrections are exactly known. 
We have demonstrated the practical effectiveness of the techniques by applying them to the II2O molecule adsorbed on the 
MgO (001) surface. We find that in this case nuclear quantum corrections to the chemical potential are only ^ 30 meV, 
but we have noted the possibility of applying the techniques to other systems where the corrections should be much larger. 
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Appendix: Test by adsorption on the rigid surface 

In Sec. |3j we presented some tests of the path-integral techniques applied to models in which carefully chosen external 
potentials act on the molecule in free space. We summarise here a different kind of test, in which full path-integral 
calculations of the mean force as a function of height z are used to calculate the chemical potential for surface adsorption, 
but with the MgO slab held rigid. These calculations are of exactly the same kind as the full calculations reported in Sec.|4j 
in that the entire system (MgO slab -I- H2O molecule) is treated by DFT. However, we modify the system to ensure that 
it is closely harmonic with the molecule on the surface by adding artificial potentials of the kind already used in Sec. [3j 
this enables us to obtain quasi-exact values for the difference between the quantum and classical values of the chemical 
potential, which can be used to validate the path-integral methods. 

The potentials that we have added to the DFT total energy to ensure almost harmonic behaviour have the form 



Kxt + llcxt + Kcyit- Here, T4xt is the potential acting on the normal to the molecular plane (Sec. 3.2), and Wcxt is the 



potential acting on the molecular bisector (Sec. 3.4). The potential /^cxt is a harmonic potential confining the motion of 



the molecular centre of mass in the plane parallel to the surface: 



ifext = + {y"'T) , (20) 



with r'^™ the centre-of-mass position defined in Sec. |3.1| The Kxt potential is zero when the normal to the molecular plane 
points along the z-axis, which in the present context means the normal to the MgO (001) surface. The Woxt potential is 
zero when the molecular bisector points along the x-axis, and we take this to be a cubic axis of the crystal parallel to the 
surface. Finally, we take the iCext potential to be zero at a point near a surface Mg ion. All three potentials Kxt, Wext 
and Kcxt are invariant under translation of the molecule along the z-direction, and they act on the molecule in exactly the 
same way when it is far from the surface as when it is on the surface. The values of the spring constants a, 7 and /3 are 
chosen to be 2.00 eV, 2.00 eV and 10.0 eV/A^ respectively. From the harmonic vibration frequencies of the molecule on 
the surface and in free space, we easily derive the chemical potential A/i^ at any temperature T, using either classical or 
quantum theory. At T = 100 K, we find that the difference (A^^)qc (quantum minus classical) has the small value 9 meV. 

As in Sec. |4j the potential of mean force calculations were performed with 8 beads, and the results were then corrected 
by thermodynamic integration with respect to number of beads, going from 8 to 16 to 32 to 64 beads. For comparison, the 
PMF calculations were also performed with classical nuclei. We find that with 8 beads the quantum-classical difference 
(A/i^)qc is —13 ± 1 meV, which differs substantially from the quasi-exact value of 9 meV. However the corrections from 
thermodynamic integration over the number of beads are large, and they do not cancel between gas and surface. In the gas 
phase, the change of free energy on going from 8 to 64 heads is 198 ±3 meV, while on the surface the change is 180 ±4 meV. 
This means that going from 8 to 64 beads stabilises the molecule on the surface by 18 meV, and the resulting corrected 
value for (A^^)qc comes to 5 ± 5 meV, which agrees well with the quasi-exact value within statistical errors. 
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